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High-accuracy binary black hole simulations are presented for black holes with spins anti-aligned 
with the orbital angular momentum. The particular case studied represents an equal-mass binary 
with spins of equal magnitude S/m^ = 0.43757± 0.00001. The system has initial orbital eccentricity 
~ 4 X 10~^, and is evolved through 10.6 orbits plus merger and ringdown. The remnant mass and 
spin are Mf = (0.961109 ± 0.000003)M and Sf/Mf^ = 0.54781 ± 0.00001, respectively, where M 
is the mass during early inspiral. The gravitational waveforms have accumulated numerical phase 
errors of < 0.1 radians without any time or phase shifts, and < 0.01 radians when the waveforms 
are aligned with suitable time and phase shifts. The waveform is extrapolated to infinity using 
a procedure accurate to < 0.01 radians in phase, and the extrapolated waveform differs by up to 
0.13 radians in phase and about one percent in amplitude from the waveform extracted at finite 
radius r — 350M. The simulations employ different choices for the constraint damping parameters 
in the wave zone; this greatly reduces the effects of junk radiation, allowing the extraction of a clean 
gravitational wave signal even very early in the simulation. 

PACS numbers: 04.25.D-, 04.25. dg, 04.30.-w, 04.30. Db, 02.70.Hm 
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I. INTRODUCTION 

Much progress has been made in recent years in the 
numerical solution of Einstein's equations for the inspi- 
ral, merger, and ringdown of binary black hole systems. 
Since the work of Pretorius [l| and the development of 
the moving puncture method 0, 0], numerical simula- 
tions have been used to analyze post-Newtonian approx- 
imations I, i, 1 0, i, i M M, M, M M M, M E M 

[Tgj . to investigate the recoil velocity of the final black 

hole J23jpj22^^ m, m, [13, M M M M M 

m, |34 l35l. l36l. l37l. iBSfT and to explore the the orbital 
dynamics of spinning binaries [H, |3^, [13, El, IH, . 

Numerical simulations can provide an accurate knowl- 
edge of gravitational waveforms, which is needed to make 
full use of the information obtained from gravitational- 
wave detectors such as LIGO and LISA. Not only can 
detected gravitational waveforms be compared with nu- 
merical results to measure astrophysical properties of 
the sources of gravitational radiation, but the detection 
probability itself can be increased via the technique of 
matched filtering (43 | , in which noisy data are convolved 
with numerical templates to enhance the signal. 

The production of accurate numerical waveforms is 
computationally expensive, making it challenging to con- 
struct an adequate waveform template bank covering a 
sufficiently large region of the parameter space of black 
hole masses and spins. One way of increasing efficiency 
is to adopt techniques known as spectral methods. For 
smooth solutions, spatial discretization errors of spectral 
methods decrease exponentially with increasing numer- 
ical resolution. In contrast, errors decrease polynomi- 
ally for the finite difference methods used in most binary 
black hole simulations. Not only have spectral methods 
been used to prepare very accurate initial data [i^, H^, 



m m m m, M m, M M, m, M 153, issi, Imi , but they 

have been used to generate the longest and most accurate 
binary black hole simulation to date [60j . 

Following the previous work of [60t , this paper presents 
the first spectral simulation of an orbiting and merging 
binary with spinning black holes: an equal mass system 
with spins of the black holes anti-aligned with the or- 
bital angular momentum. Simulations of binaries with 
spins parallel to the orbital momentum are certainly not 
new, e.g. [ill, HI, H [H, [H, lH, [6l|. Our goal here is 
to show that such systems can be simulated with spec- 
tral methods, and that the high accuracies achieved for 
the non-spinning case carry over into this more general 
regime. 

The spin of each black hole is S/m'^ = 0.43757 ± 
0.00001. The determination of this quantity, as well 
as other spin measures, is explained in more detail in 
Sec. IIVBI The evolution consists of 10.6 orbits of in- 
spiral with an orbital eccentricity of e 4 x 10~^, fol- 
lowed by the merger and ringdown. We find that this 
simulation has accuracy comparable to that of the simu- 
lation presented in [g^I ■ We also present different choices 
for the constraint damping parameters in the wave zone; 
these choices cause the initial noise ("junk radiation") to 
damp more rapidly, resulting in a useable, almost noise- 
free waveform much earlier in the simulation. 

This paper is organized as follows: In Sec. HIl we dis- 
cuss the construction of our initial data. In Sec. IIIIl we 
describe the equations, gauge conditions, and numerical 
methods used to solve Einstein's equations. In Sec. IIVI 
we present several properties of our simulations, includ- 
ing constraints, and the spins and masses of the black 
holes. In Sec. |Vl we explain the extraction of gravita- 
tional waveforms from the simulation, and the extrapo- 
lation of the waveforms to infinity. Finally, in Scc. lVIi we 



2 



discuss outstanding difficulties and directions for future 
work. 



II. INITIAL DATA 

The initial data are almost identical to those used 
in the simulation of an equal-mass, non-spinning black 
hole binary presented in Refs. M. IGOI . We use quasi- 
equilibrium initial data [s^, Issl. |62|| (see also [isl. l49|') 
built using the conformal thin sandwich formalism |63l . 
H^l, and employing the simplifying choices of conformal 
flatness and maximal slicing. Quasi-equilibrium bound- 
ary conditions arc imposed on spherical excision bound- 
aries for each black hole, with the lapse boundary condi- 
tion given by Eq. (33a) of Ref. 
are centered at 

and C2 = (— c?/2, 0, 0), where we choose the same coordi 
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The excision spheres 
Cartesian coordinates CI = (d/2,0,0) 



nate distance d and the same excision radii as in 

Within this formalism, the spin of each black hole is 
determined by a parameter fir and a conformal Killing 
vector (tangential to the excision sphere); these enter 
into the boundary condition for the shift at an excision 
surface 1531 . We will use the sign convention of Eq. (40) 
in Ref. [57j , so that positive f2r corresponds to corotating 
black holes. The same value of ilr is chosen at both ex- 
cision surfaces, resulting in black holes with equal spins. 
In Refs. [1, H^, was chosen to ensure that the black 
hole spins vanish (S^- In this paper, we instead fix fir 
at some negative value, resulting in moderately spinning 
black holes that countcrrotate with the orbital motion. 

Two more parameters need to be chosen before initial 
data can be constructed: The orbital angular frequency 
flo and the radial velocity Vr of each black hole. These 
parameters are determined by an iterative procedure that 
minimizes the orbital eccentricity during the subsequent 
evolution of the binary: We start by setting flo and Vr 
to their values in the non-spinning evolution of Ref. [6^ , 
we solve the initial value equations a pseudo-spectral el- 
liptic solver [5l|, and we evolve for about 1-2 orbits us- 
ing the techniques described in Sec. IIIII Analysis of this 
short evolution yields an estimate for the orbital eccen- 
tricity, and improved parameters flo and Vr that result in 
a smaller orbital eccentricity. This procedure is identical 
to Ref. 3, except that we include a term quadratic in t 
for the function used to fit the radial velocity (ds/dt), 
to obtain better fits. We repeat this procedure until 
the eccentricity of the black hole binary is reduced to 
e ~ 4 X 10~^. Properties of this low-ccccntricity initial 
data set are summarized in the top portion of Tabic [H 

The data in the upper part of Table [J arc given in 
units of Mid, the sum of the black hole masses in the 
initial data. For any black hole (initial data, during the 
evolution, the remnant black hole after merger), we define 
its mass using Christodoulou's formula. 



Initial data 




Coordinate separation 


d/MiD = 13.354418 


Radius of excision spheres 


rexc/A/iD = 0.382604 


Orbital frequency 


noMiD = 0.0187862 


Radial velocity 


Vr = -7.4710123 X 10"^ 


Orbital frequency of horizons 


QrMm = -0.242296 


Black hole spins 


XiD = 0.43785 


ADM energy 


Madm/Mid = 0.992351 


Total angular momentum 


JadmM^d = 0.86501 


Initial proper separation 


So /A/id = 16.408569 


Evolution 




Initial orbital eccentricity 


e ?a 4 X lO"-'^ 


Mass after relaxation 


M = (1.000273 ± 0.000001)A/iD 


Spins after relaxation 


X = 0.43757 ± 0.00001 


Time of merger (common AH) tcAu = 2399.38A'f 


Final mass 


Mf = (0.961109 ± 0.000003)A/ 


Final spin 


Xf = 0.54781 ± 0.00001 



TABLE I: Summary of the simulation presented in this paper. 
The first block lists properties of the initial data, the second 
block lists properties of the evolution. 



We use the apparent horizon area ^ah to define the irre- 
ducible mass TOirr = \/ ^AH/(167r). The nonnegative spin 
S of each black hole is computed with the spin diagnos- 
tics described in [s^ . Unless noted otherwise, we com- 
pute the spin from an angular momentum surface inte- 
gral [66l . [U using approximate Killing vectors of the ap- 
parent horizons, as described in [KJ,!!!] (see also [69l.[70l|). 
We define the dimcnsionlcss spin by 



S 



(2) 



m 



4mi, 



(1) 



III. EVOLUTIONS 

A. Overview 

The Einstein evolution equations are solved with the 
pseudo-spectral evolution code described in Ref. 60| . 
This code evolves a first-order representation [71 | of 
the generalized harmonic system [tI, [tI, [tII and in- 
cludes terms that damp away small constraint viola- 
tions [71I, [zl, [ill. The computational domain extends 
from excision boundaries located just inside each appar- 
ent horizon to some large radius, and is divided into sub- 
domains with simple shapes (e.g. spherical shells, cubes, 
cylinders). No boundary conditions are needed or im- 
posed at the excision boundaries, because all character- 
istic fields of the system are outgoing (into the black hole) 
there. The boundary conditions on the outer bound- 
ary [71I, [7^, [73 are designed to prevent the influx of un- 
physical constraint violations [zl. 111, HO, [Ml, [83. l83l . l83 
and undesired incoming gravitational radiation [8 a . l86l | , 
while allowing the outgoing gravitational radiation to 
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Run 


Npts CPU-h CPU-h/T 


Nl 


(64^65^65^) 9,930 3.4 


N2 


(70^72^72^) 16,195 5.6 


N3 


(76■^78^80=*) 28,017 9.7 


N4 


(82^84^873) 44,954 15.5 



TABLE II: Approximate number of collocation points and 
CPU usage for the evolutions performed here. The first col- 
umn indicates the name of the run. Npts is the approximate 
number of collocation points used to cover the entire com- 
putational domain. The three values for Npts are those for 
the inspiral, plunge, and ringdown portions of the simula- 
tion, which are described in Sections IIII CI IIII Dl and IIII El 
respectively. The total run times T are in units of the total 
Christodoulou mass M [cf. Eq. ([T}] of the binary. 



pass freely through the boundary. Interdomain boundary 
conditions are enforced with a penalty method [13, Ull . 

The gauge freedom in the generalized harmonic formu- 
lation of Einstein's equations is fixed via a freely speci- 
fiable gauge source function Ha that satisfies the con- 
straint 



— Ca — Tab'' + Ha 



(3) 



where V^bc are the spacetime Christoffel symbols. We 
choose Ha differently during the inspiral, plunge and 
ringdown, as described in detail in Sections IIII Ci IIII Dl 
and [mil 

In order to treat moving holes using a fixed grid, 
we employ multiple coordinate frames [89|: The equa- 
tions are solved in 'inertial frame' that is asymptotically 
Minkowski, but the grid is fixed in a 'comoving frame' 
in which the black holes do not move. The motion of 
the holes is accounted for by dynamically adjusting the 
coordinate mapping between the two frames^. This co- 
ordinate mapping is chosen differently at different stages 
of the evolution, as described in Sections IIII CI IIII Dl 
and [mil 

The simulations are performed at four different resolu- 
tions, Nl to N4. The approximate number of collocation 
points for these resolutions is given in Table [III 



B. Relaxation of Initial Data 

The initial data do not precisely correspond to two 
black holes in equilibrium, e.g., because tidal deforma- 
tions are not incorporated correctly, and because of the 
simplifying choice of conformal flatness. Therefore, early 
in the evolution the system relaxes and settles down into 
a new steady-state configuration. Figure [1] shows the 
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^ All coordinate quantities (e.g. trajectories, waveform extraction 
radii) in this paper are given with respect to the inertial frame 
unless noted otherwise. 
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FIG. 1: Irreducible mass (top panel) and spin (bottom panel) 
of the black holes during the relaxation of the initial data to 
the equilibriun (steady-state) inspiral configuration. Shown 
are four different numerical resolutions, Nl (lowest) to N4 
(highest), cf. Table [iTl Up to t ~ lOM, both mass and spin 
change by a few parts in 10*, then they remain approximately 
constant (as indicated by the dashed horizontal lines) until 
shortly before merger. These steady-state values are used to 
define M and X- 



change in irreducible mass and spin relative to the ini- 
tial data during the evolution. During the first ~ lOM 
of the evolution, Afin- increases by about 3 parts in lO'' 
while the spin decreases by about 1 part in 10''. These 
changes are resolved by all four numerical resolutions, 
labeled Nl (lowest) to N4 (highest), and converge with 
increasing resolution. After the initial relaxation, for 
lOM ^t< 2350 Af, the mass is constant to about 1 part 
in 10^, as can be seen from the convergence of the differ- 
ent resolutions in the upper panel of Fig. [1] In the last 
~ 5GA/ before merger, the mass increases slightly (seen as 
a vertical feature at the right edge of the plot), an effect 
we will discuss in more detail in the context of Fig.[5l The 
spin is likewise almost constant for 10 < t/M < 1000, al- 
though some noise is visible for t < lOOA/. 

We shall take the steady-state masses and spins eval- 
uated at t ~ 200Af as the physical parameters of the 
binary being studied. Specifically, all dimcnsionful c^uan- 
titics will henceforth be expressed in terms of the mass 
scale A/, which we define as the total mass after relax- 
ation. 

The relaxation of the black holes in the first ~ lOAf 
of the evolution is also accompanied by the emission 
of a pulse of unphysical "junk radiation." This pulse 
passes through the computational domain, and leaves 
through the outer boundary after one light crossing time. 
The junk radiation contains short wavelength features, 
which are not resolved in the wave zone. It turns out 
that the constraint damping parameters 70 and 72 (see 
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[71| ) influence how the unresolved junk radiation inter- 
acts with the numerical grid. Large constraint damping 
parameters enhance the conversion of the outgoing junk 
radiation (at the truncation error level) into incoming 
modes. This incoming radiation then lingers for several 
light-crossing times within the computational domain, 
imprinting noise into the extracted gravitational radia- 
tion. For small constraint damping parameters, this con- 
version is greatly suppressed, and numerical noise due 
to junk radiation diminishes much more rapidly. The 
simulations presented here use 7o = 72 0.00225/M 
in the wave zone; these values are smaller by a fac- 
tor 100 than those used in I, iO]. (Even smaller con- 
straint damping parameters fail to suppress constraint 
violations. Note that constraint damping parameters are 
much larger, 70 = 72 '^-^ 3.56/A/. in the vicinity of the 
black holes.) The waveforms presented here show conse- 
quently reduced contaminations in the early part of the 
evolution will be discussed in Sec. IVll cf. Fig. [T5l 



C. Inspiral 

During the inspiral, the mapping between the comov- 
ing and inertial frames is chosen in the same way as in 
Refs. HOl and is denoted by A4i : x'* a;*, where 
primed coordinates denote the comoving frame and un- 
primed coordinates denote the inertial frame. Explicitly, 
this map is 



a{t) + {1 ~ a{t)) 

^0 . 



= 0' + bit), 



(4) 

(5) 
(6) 



where (r, 9, (f)) and (r', 9' ,(/)') denote spherical polar coor- 
dinates relative to the center of mass of the system in in- 
ertial and comoving coordinates, respectively. We choose 
i?Q = 467M. The functions a{t) and b{t) are determined 
by a dynamical control system as described in Ref. [soj . 
This control system adjusts a{t) and b{t) so that the cen- 
ters of the apparent horizons remain stationary in the 
comoving frame. 

While each hole is roughly in equilibrium during inspi- 
ral, we choose the gauge source function Ha in the same 
way as in Refs. 0,[63l: A new quantity Ha is defined that 
has the following two properties: 1) Ha transforms like 
a tensor, and 2) in inertial coordinates Ha ~ Ha- Ha is 
chosen so that the constraint Eq. Q is satisfied initially, 
and Ha' is kept constant in the comoving frame, i.e., 

dt'Ha' = 0. (7) 

Here primes refer to comoving frame coordinates. This 
is essentially an equilibrium condition. 



D. Plunge 

We make two key modifications to our algorithm to 
allow evolution through merger. The first is a change in 
gauge conditions, as in Ref. [60|. The second is a change 
in coordinate mappings that allows the excision bound- 
aries to more closely track the horizons. We describe 
both of these changes here. 

Following Ref. [60| , at some time t = tg (where g stands 
for "gauge") we promote the gauge source function Ha 
to an independent dynamical field that satisfies 



W,Ha - Qa{x, i, i^ab) + ^2t'dbHa. 



(8) 



Here V^Vc is the curved space scalar wave operator (i.e. 
each component of Ha is evolved as a scalar), ipab is the 
spacetime metric, and is the timelikc unit normal to 
the hypersurface. The driving functions Qa are 



Qt = /(x,t)ei 

Q, = g{x,t)£,3 



l-N 



(9) 
(10) 



where N and /?' are the lapse function and the shift 
vector, 77, ^1, ^2, and ^3 are constants, and f{x,t), and 
g{x, t) are prescribed functions of the spacetime coordi- 
nates. Eq. dSt is evolved in first-order form, as described 
in Ref. [60|. Eq. ^ requires values of Ha and its time 
derivative as initial data; these are chosen so that Ha and 
dfHa are continuous a,t t = tg. 

This gauge is identical to the one used in Ref. [60t . 
except that the parameters and functions that go into 
Eq. dH]) are chosen slightly differently: We set 77 = 4, 
?! = 0.1, 6 = 6.5, ^3 = 0.01, and 



X (1 V-i)e-'-'V-l 



(11) 



g{x,t) = (l_e"(*-*«'/"^) 

X (l-e-(*-*«)'/'^')(t-tg)e-'-"/'^3, (12) 

where r' is the coordinate radius in comoving coordi- 
nates, and the constants arc ai ^ 62A/, (T2 ^ 44. 5M, 
CT3 ~ 35M, (74 ~ 4.5M, and as ~ 3M. The function 
g{x,t) in Qi, which drives the shift towards zero near 
the black holes, has a factor {t — tg) that is absent in 
Ref. [131. Prescribing g{x,t) in this way drives the shift 
towards zero more strongly at late times, which for this 
case is more effective in preventing gauge singularities 
from developing. 

The second change we make a.t t = tg is to control the 
shape of each excision boundary so that it matches the 
shape of the corresponding apparent horizon. In the co- 
moving frame, where the excision boundaries are spheri- 
cal by construction, this means adjusting the coordinate 
mapping between the two frames such that the apparent 
horizons are also spherical. Without this "shape con- 
trol" , the horizons become sufficiently distorted with re- 
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FIG. 2: Coordinate trajectories of the centers of the apparent 
horizons represented by the blue and red curves, up until the 
formation of a common horizon. The closed curves show the 
coordinate shapes of the corresponding apparent horizons. 

spect to the excision boundaries that the excision bound- 
aries fail to remain outflow surfaces and our excision al- 
gorithm fails. For the non-spinning black hole binary in 
Ref. (60j . shape control was not necessary before merger. 
To control the shape of black hole 1, we define the map 
Mari ■■ x" x\ 

6 ^ 0', (13) 
- 0', (14) 

f ^ r'-q,ir')J2 E >^Lme..{0' ,<P'), (15) 

1=0 m=-l 

where 

gi(r') =e-('''-'-«(*»'/"', (16) 

and {r',9',(l)') are spherical polar coordinates centered 
at the (fixed) comoving-coordinatc location of black hole 
1. The function qi{r') limits the action of the map to 
the vicinity of hole 1. The constant cr^ is chosen to be 
r-^ A.5M, and r'^it) = rg-t-z^i(t — tg)^-^ is a function of time 
that approximately follows the radius of the black hole, 
with constants Tq ~ 1.2M and i^i ~ 0.00046M. Similarly, 
we define the map A1ah2 for black hole 2. Then the full 
map A^ni : — > from the comoving coordinates x"' 
to the inertial coordinates is given by 

Mm-^ MioMah2oMahi- (17) 

The functions Xl^{t) and Xj^^it) arc determine d by dy- 
namical control systems as described in Refs. [60, [89| . 
so that the apparent horizons are driven to spheres (up 



to spherical harmonic component I = ^max) in comoving 
coordinates. Note that A^ah i : x'* is essentially 

the same map that we use to control the shape of the 
merged horizon during ringdown, and the control system 
for that map (and for the map A^ah2) is the same as the 
one described in Ref. [gO] for controlling the shape of the 
merged horizon. 

In addition to the modifications to the gauge condi- 
tions and coordinate map described above, the numer- 
ical resolution is also increased slightly around the two 
black holes during this more dynamical phase, and the 
evolution is continued until time tm, shortly after the 
formation of a common horizon. The coordinate trajec- 
tories of the apparent horizon centers are shown in Fig. [2] 
up until tm, at which point the binary has gone through 
10.6 orbits. 



E. Ringdown 

Our methods for continuing the evolution once a com- 
mon horizon has formed are the same as in Ref. (60j . Af- 
ter a common apparent horizon is found, all variables are 
interpolated onto a new computational domain that has 
only a single excised region. Then, a new comoving coor- 
dinate system (and a corresponding mapping to inertial 
coordinates) is chosen so that the new excision boundary 
tracks the shape of the apparent horizon in the inertial 
frame, and also ensures that the outer boundary behaves 
smoothly in time. The gauge conditions are modified as 
well: the shift vector is no longer driven to zero, so that 
the solution can relax to a time-independent state. This 
is done by allowing the gauge function g{x, t) that ap- 
pears in Eq. (fTUl) to gradually approach zero; the gauge 
source function Ha still obeys Eqs. ([Sl fTU]) as during the 
plunge. Specifically, we change the functions f{x,t) and 
g[x, t) from Eqs. dTT]) and ^ to 

/(x,t) = (2-e-(*-*«)/"0 

X (l-e"(*-*«)'/"')e-'^"'/"3, (18) 

g{x,t) = (1 -e-^*-*')/'^^) 

X (1 - e-(*-*«)'/'^')(t - tg)e~''"'/'"3 

X e-(*-*'")/'"«, (19) 

where r" is the coordinate radius in the new comov- 
ing coordinates, erg ~ S.lAf, and (here m stands for 
"merger") is the time we transition to the new domain 
decomposition. 

IV. PROPERTIES OF THE NUMERICAL 
SOLUTIONS 

A. Constraints 

We do not explicitly enforce either the Einstein con- 
straints or the secondary constraints that arise from writ- 
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FIG. 3: Constraint violations of runs on different resolutions. 
The top panel shows the norm of all constraints, normal- 
ized by the norm of the spatial gradients of all dynamical 
fields. The bottom panel shows the same data, but without 
the normalization factor. The norms are taken over the 
portion of the computational volume that lies outside appar- 
ent horizons. Note that the time when we change the gauge 
before merger, tg ~ 2370Af , and the time when we regrid 
onto a new single-hole domain after merger, tm ~ 2400Af , 
are slightly different for different resolutions. 
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FIG. 4: Dimensionless spins x of one black hole in the N4 
evolution, evaluated using an approximate Killing vector, a 
coordinate rotation vector — 9^, or the extrema of the instrin- 
sic scalar curvature on the apparent horizon. Bottom panels 
show detail at early and late times. Also shown are the time 
of gauge change tg before merger, and the time tm that we 
transition to a single-hole evolution just after merger. 



cnt horizons, and this region is newly excised from the 
computational domain at t = tm- 



ing the system in first-order form. Therefore, examining 
how well these constraints are satisfied provides a useful 
consistency check. Figure |3] shows the contraint viola- 
tions for the evolutions at different resolutions. The top 
panel shows the norm of all the constraint fields of 
our first-order generalized harmonic system, normalized 
by the norm of the spatial gradients of the dynamical 
fields (see Eq. (71) of Ref.[2l|)- The bottom panel shows 
the same quantity, but without the normalization factor 
(i.e., just the numerator of Eq. (71) of Ref.[ll|)- The 
norms are taken over the portion of the computational 
volume that lies outside the apparent horizons. 

The constraints increase as the black holes approach 
each other and become increasingly distorted. At tg = 
2372.05M for N4 {tg = 2372.05M for N3, tg = 2376.5M 
for N2, tg = 2376. 5M), the gauge conditions are changed 
(cf. 1111 Pp and the resolution around the holes is in- 
creased slightly. Because of the change in resolution, 
the constraints drop by more than an order of magni- 
tude. Close to merger, the constraints grow larger again. 
The transition to a single-hole evolution (cf. 1111 Ep oc- 
curs at t,n = 2399.647\/ for N4 {tm = 2399.66M for N3, 
tm = 2401.27M for N2, tm = 2404.23M for Nl). Shortly 
after this time, the constraints drop by about two or- 
ders of magnitude. This is because the largest constraint 
violations occur near and between the individual appar- 



B. Black hole spins and masses 

There are different ways to compute the spin x(i) of 
a black hole. The approach we prefer computes the spin 
from an angular momentum surface integral [66l . [6^ using 
approximate Killing vectors of the apparent horizons, as 
described in [s^, Hal (see also [bI, [TO] ) . We shall denote 
the resulting spin by xakv(0- Another less sophisticated 
method simply uses coordinate rotation vectors, and we 
denote the resulting spin by XCoord(0- We also use two 
more spin diagnostics that are based on the minimium 
and maximum of the instrinsic scalar curvature of the 
apparent horizon for a Kerr black hole [s^l ; we call these 
Xsc"(0 and Xsc'^(^)- These last two measures of spin 
are expected to give reasonable results when the black 
holes are sufficiently far apart and close to equilibrium, 
and after the final black hole has settled down to a time- 
independent state. However, they are expected to be less 
accurate near merger and at the start of the evolution. 

Fig. |4] shows these four spin measures for black hole 
1 in the N4 evolution during inspiral and plunge. From 
the lower left panel we see that Xsc"(^) ''^^'•^ Xsc'^W 
fer from XCoord(i) and XAKv(i) by more than a factor of 
two a.t t ~ 0. This indicates that the initial black holes 
do not have the appropriate shape for the Kerr solution; 
i.e. they are distorted because of the way the initial data 
is constructed. As the black holes relax, Xsc"(^) ^^'^ 
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FIG. 5: Sum of Christodoulou masses M{t) and sum of irre- 
ducible masses Mirr(t) of the two black holes during inspiral. 
The data is from the N4 evolution, and uses Xakv when com- 
puting M{t). Insets show detail at late times, and indicate 



the transition times tg and t„ 



Xsc^i^) approach the other two spin measures. The re- 
laxed spin at t - 200Af is x = 0.43757 ± 0.00001, where 
the uncertainty is based on the variation in xakv be- 
tween t = lOOM and t = lOOOM. During the inspiral, 
XAKv{t) decreases slowly and monotonically, dropping 
by 10~* at 90M before merger, and dropping by 0.01 at 
the time of merger. Tidal dissipation should slow down 
the black holes, so this decrease is physically sensible. In 
contrast, the other three spin diagnostics show a mild 
increase in spin, suggesting that they arc less reliable. 
Close to merger, Xsc"('^) and Xsc'^('') increase dramat- 
ically, with Xsc'^('^) growing as large as 0.92. In this 
regime, the shapes of the individual black holes are dom- 
inated by tidal distortion, and are therefore useless for 
measuring the spin. 

The Christodoulou mass m of one black hole, as de- 
fined in Eq. ([1]), depends on the spin. We take xakv(^) 
as the preferred spin measure, and use it to compute the 
total Christodoulou mass M{t) during the inspiral and 
plunge. This is shown in the top panel of Fig. [S] The 
Christodoulou mass settles down to M{t)/M = 1.000000 
after t = 150M (this defines M), and increases to 
M{t)/M = 1.00114 at the time of merger. Most of the 
increase in mass occurs very close to merger, as can be 
seen from the inset of Fig. [S] Until about 30M before 
merger (i.e. t = 2370M), the mass is constant to a few 
parts in 10®. For comparison, in the bottom panel we also 
display Mi„{t), the sum of the irreducible masses, which 
does not depend on the spin. This quantity settles down 
to Mi„{t)/M = 0.974508 at t = 200M, and increases to 
Mi„{t)/M = 0.97668 at t = 2400Af. Again, almost all of 
this increase happens shortly before merger. During the 
inspiral up to 30M before merger, Mi„{t)/M increases 
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FIG. 6: Dimensionless spins x(^) of the final black hole 
in the N4 evolution. The (most reliable) spin diagnostic 
Xakv starts at ~ 0.4 and increases to its final value Xf = 
0.54781 ± 0.00001. The other spin diagnostics are unreliable 
for the highly distorted black hole shortly after merger, but 
subsequently approach xakv. 



by only 6 x 10 ^, but in the last 30A/ the increase is 
- 0.002. 

The merger results in one highly distorted black hole, 
which subsequently rings down into a stationary Kerr 
black hole. Figure [5] shows our four spin diagnostics dur- 
ing the ringdown. The spin measures Xsc"(0 and Xsc''('') 
assume a Kerr black hole. Just after merger, the horizon 
is highly distorted, so these two spin diagnostics are not 
valid there. However, as the remnant black hole rings 
down to Kerr, Xsc'^(^) Xsc"(^) approach the quasi- 
local AKV spin to better than 1 part in 10^ (see the inset 
of Fig. [6|). The quasi- local spin based on coordinate ro- 
tation vectors, XCoord(Oi also agrees with the other spin 
measures to a similar level at late times. The spin of 
the final black hole points in the direction of the initial 
orbital angular momentum. 

The Christodoulou mass Mf{t) of the final black hole 
in the N4 evolution, again evaluated using xakv(^)i is 
shown in the top panel of Fig. [T] The mass settles down 
to a final value of Mf/M = 0.961109 ± 0.000003. The 
bottom panel shows the irreducible mass Mh-r j(<) of the 
final black hole, which settles down to a final value of 
Mirr.f = 0.921012 ± 0.000003. The uncertainties are de- 
termined from the difference between runs N4 and N3, 
so they include only numerical truncation error and not 
any systematic effects. The uncertainty in the mass is 
visible in the insets of Fig. [T] 
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FIG. 7; The top panel shows the Christodoulou mass Mf{t) 
of the final black hole in the N4 and N3 runs, computed us- 
ing XAKv(t). The bottom panel shows the irreducible mass 
Mi„,f(t). 



V. COMPUTATION OF THE WAVEFORM 

A. Waveform extraction 

Gravitational waves are extracted from the simulution 
on spheres of different values of the coordinate radius 
r, following the same procedure as in Refs. f60l. [65l. loot . 
The Newman- Penrose scalar ^'4 in terms of spin- weighted 
spherical harmonics of weight -2: 

VI/4 {t, r, 9,ct,)^Y. *4™ (t. 0-2 Yim. {0, </)) , (20) 

Im 

where the ^'4™ arc expansion coefficients defined by this 
equation. Here we also focus on the dominant (l,m) = 
(2, 2) mode, and split the extracted waveform into real 
phase ^ and real amphtude A, defined by (see e.g. [l,[9l|) 

*f (r,t) = A(r,t)e-"^('"^*). (21) 

The gravitational-wave frequency is given by 



dt 



(22) 



The minus sign in Eq. (|2ip is chosen so that the phase 
increases in time and to is positive. 

The coordinate radius of our outer boundary is located 
at i?i„ax = 427M at t ^ and i?,„ax = 365M at t > 
2500M; it shrinks slightly during the evolution because 
of the mappings [cf. Eq. H])] used in our dual frame 
approach. The (/, m) = (2, 2) waveform, extracted at a 
single coordinate radius r = 350i\/ for the N4 evolution, 
is shown in Fig [51 The short pulse at t ~ 360M is due 
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FIG. 8: Gravitational waveform extracted at finite radius r = 
350M for the N4 evolution. The top panel zooms in on the 
inspiral waveform, and the bottom panel zooms in on the 
merger and ringdown. 



to junk radiation. The magnitude of this pulse is about 
twice as large as for non-spinning black holes, cf. Ref. [1, 



B. Convergence of extracted waveforms 

In this section we examine the convergence of the grav- 
itational waveforms extracted at fixed radius, without 
extrapolation to infinity. This allows us to study the 
behavior of our code without the complications of ex- 
trapolation. The extrapolation process and the resulting 
extrapolated waveforms are discussed in Sec. IV CI 

Figure [9] shows the convergence of the gravitational- 
wave phase 4) and amplitude A with numerical resolu- 
tion. For this plot, the waveform was extracted at a 
fixed inertial-coordinate radius of r = 350M. Each line 
in the top panel shows the absolute difference between 
(j) computed at some particular resolution and </> com- 
puted from our highest resolution N4 run. The curves 
in the bottom panel similarly show the relative ampli- 
tude differences. When subtracting results at different 
resolutions, no time or phase adjustment has been per- 
formed. The noise at early times is due to junk radiation 
generated near t = 0. Most of this junk radiation leaves 
through the outer boundary after one crossing time. The 
plots show that the phase difference accumulated over 
10.6 orbits plus merger and ringdown — in total 31 grav- 
itational wave cycles — is less than 0.1 radians, and the 
relative amplitude differences are less than 0.017. These 
numbers can be taken as an estimate of the numerical 
truncation error of our N3 run. Because of the rapid 
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FIG. 9: Convergence of gravitational waveforms with nu- 
merical resolution. Shown are phase and amplitude differ- 
ences between numerical waveforms ^'4^ computed using dif- 
ferent numerical resolutions. All waveforms are extracted at 
r = 350M, and no time shifting or phase shifting is done to 
align waveforms. 




FIG. 10: Convergence of gravitational waveforms with numer- 
ical resolution. Same as Fig.|9]except all other waveforms are 
time-shifted and phase-shifted to best match the waveform of 
the N4 run. 



convergence of the code, we expect that the errors of the 
N4 run are significantly smaller. 

Figure [To] is the same as Fig. [9] after the Nl, N2, N3 
waveforms have been time shifted and phase shifted to 
best match the waveform of the N4 evolution. This best 
match is determined by a simple least-squares procedure: 



we minimize the function 
J2 (v4l(^^)e**l(*') - A2(i.-f ^o)e*^*'^*'+*°^+*"^)^ (23) 

i 

by varying to and (j)o ■ Here Ai , 0i , A2 , and 02 are the am- 
plitudes and phases of the two waveforms being matched, 
and the sum goes over all times ti at which waveform 
1 is sampled. This type of comparison is relevant for 
analysis of data from gravitational-wave detectors: when 
comparing experimental data with numerical detection 
templates, the template will be shifted in both time and 
phase to best match the data. For this type of compar- 
ison, Fig. [TU] shows that the numerical truncation error 
of our N3 run is less than 0.01 radians in phase and 0.1 
percent in amplitude for t > 550A/. At earlier times, the 
errors are somewhat larger and are dominated by residual 
junk radiation. 



C. Extrapolation of waveforms to infinity 

Gravitational-wave detectors measure waveforms as 
seen by an observer effectively infinitely far from the 
source. Since our numerical simulations cover only a 
finite spacetime volume, after extracting waveforms at 
multiple finite radii, we extrapolate these waveforms to 
infinite radius using the procedure described in [60j (see 
also [lOl for more details). This is intended to reduce 
near-field effects as well as gauge effects that can be 
caused by the time dependence of the lapse function or 
the nonoptimal choice of tetrad for computing 4*4. 

The extrapolation of the extracted waveforms involves 
first computing each extracted waveform as a function 
of retarded time u = ts — r* and extraction radius rarcai 
(see [gOI for precise definitions). Then at each value of 
w, the phase and amplitude are fitted to polynomials in 

1/^arcal- 

<^(u,ra.eal) = -^(O) (") + E , (24) 

^arcal 

ra.cai = A^o){n) + j2^¥^- (25) 

^arcal 

The phase and amplitude of the desired asymptotic wave- 
form arc thus given by the leading-order term of the cor- 
responding polynomial, as a function of retarded time: 

(/)(u) = 0(0) (w), (26) 

rarcai A(u) = ^(0)(u)- (27) 

Figure [Tl] shows phase and amplitude differences be- 
tween extrapolated waveforms that are computed us- 
ing different values of polynomial order n in Eqs. (|24[) 
and ()25p . The extrapolation is based on waveforms ex- 
tracted at 20 different radii between 75M and 350 A'/. As 
in [60j . our preferred extrapolation order is n = 3, which 
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FIG. 11: Convergence of extrapolation to infinity for extrapo- 
lation of order n. For each n, plotted is the extrapolated wave- 
form from N4 using order n + 1 minus the extrapolated wave- 
form using order n. The top panel shows phase differences, 
the bottom panel shows amplitude differences. No shifting in 
time or phase has been done for this comparison. 
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FIG. 12: Late-time phase convergence of extrapolation to in- 
finity. Same as the top panel of Fig. 1111 except zoomed to 
late times. The peak amplitude of the waveform occurs at 
ts-r* = 2410.6M. 
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FIG. 13: Phase and relative amplitude differences between 
extrapolated and extracted waveforms for N4. The extracted 
waveform is extracted at coordinate radius r = 350M. The 
waveforms are time-shifted and phase-shifted to produce the 
best least-squares match. 



between our preferred extrapolated waveform using n = 3 
and the waveform extrapolated at coordinate radius r = 
350M, both for the N4 run. The extrapolated waveform 
has been shifted in time and phase so as to best match the 
71 = 3 extrapolated waveform, using the least-squares fit 
of Eq. ([23|) . The phase difference between extrapolated 
waveform and waveform extracted at r = 350M becomes 
as large as 0.13 radians, and the amplitude difference is 
on the order of 1 per cent. 

Figure [14] presents the final waveform after extrapola- 
tion to infinite radius. There are 22 gravitational-wave 
cycles before the maximum of |^'4|, and 9 gravitational- 
wave cycles during ringdown, over which the amplitude 
of 1 1 drops by four orders of magnitude. 



VI. DISCUSSION 



gives a phase error of less than 0.004 radians and a rel- 
ative amplitude error of less than 0.006 during most of 
the inspiral, and a phase error of less than 0.01 radians 
and a relative amplitude error of 0.006 in the ringdown. 

Figure [T^ is the same as the top panel of Fig. [Til ex- 
cept zoomed to late times. During merger and ringdown, 
the extrapolation procedure does not converge with in- 
creasing extrapolation order n: the phase differences are 
slightly larger for larger n. This was also seen for the 
extrapolated waveforms of our equal-mass nonspinning 
black hole binary and is possibly due to gauge ef- 
fects that do not obey the fitting formulae, Eqs. ((M]) 
and ((25|) . 

Figure [T^] shows the phase and amplitude differences 



We have presented the first spectral computation of 
a binary black hole inspiral, merger, and ringdown with 
spinning black holes, and find that we can achieve sim- 
ilar accuracy for the final mass, final spin, and gravita- 
tional waveforms as in the non-spinning case [60|. For 
initial spins of x = 0.43757 ± 0.00001, the mass and spin 
of the final hole are Mj/M = 0.961109 ± 0.000003 and 
Xf = 0.54781 ± 0.00001. The uncertainties are based on 
comparing runs at our highest two resolutions, and do 
not take into account systematic errors (e.g. the pres- 
ence of a finite outer boundary or gauge effects). Note 
that for the non-spinning case [60| , we found that chang- 
ing the outer boundary location produced a smaller ef- 
fect on the final mass and spin than changing the reso- 
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lution, and that the outer boundary for the evolutions 
presented here is more distant (at late times, when most 
of the radiation passes through the boundary) than it 
was in Ref. [60| . The uncertainties in the gravitational 
waveforms are < 0.01 radians in phase and < 0.6 per- 
cent in amplitude (when waveforms are time and phase 
shifted). These uncertainties are based on comparisons 
between our two highest resolution runs and comparisons 
between different methods of extrapolating waveforms to 
infinite extraction radius. 

The methods used here to simulate plunge and ring- 
down are similar to those in Ref. [60| . The primary disad- 
vantage of these methods is that they require fine tuning 
during the plunge (Sec. IIIIPp . For example, the func- 
tion g{x,t) defined in Eq. p2|) must be chosen carefully 
or else the simulation fails shortly (a few M) before a 
common horizon forms. There are at least two reasons 
that fine tuning is currently necessary. First, the gauge 
conditions must be chosen so that no coordinate singu- 
larities occur before merger. Second, the excision bound- 
aries do not coincide with the apparent horizons, but 
instead they lie somewhat inside the horizons. If the ex- 
cision boundaries exactly followed the horizons, then the 
characteristic fields of the system would be guaranteed 
to be outflowing (into the holes) at the excision bound- 
aries, so that no boundary condition is required there. 
But for excision boundaries inside the horizons, the out- 
flow condition depends on the location of the excision 
boundary, its motion with respect to the horizon, and 
the gauge. Indeed, the most common mode of failure 



for improperly-tuned gauge parameters is that the out- 
flow condition fails at some point on one of the excision 
boundaries. We have been working on improved gauge 
conditions [92| and on improved algorithms for allowing 
the excision boundary to more closely track the apparent 
horizon. These and other improvements greatly reduce 
the amount of necessary fine tuning and allow mergers 
in generic configurations, and will be described in detail 
elsewhere f93| . 

Another quite important improvement lies in the 
choice of constraint damping parameters. To illustrate 
this effect, Fig. [T5] compares the gravitational wave phase 
extrapolation for the simulation presented here with the 
similar plot for an earlier run Q with different constraint 
damping parameters. As can be seen in Fig. [121 the 
improved constraint damping parameters result in sig- 
nificantly reduced noise. For the earlier simulation, the 
waveform was unusable for t — r*< lOOOAf , and was still 
noticeably noisy at lOOOA/ < t - r* < 2000M. For the 
new simulation, the smaller constraint damping parame- 
ters result in clean waveforms as early as i — r* ~ 250Af , 
despite the observation that the spinning black holes 
result in a pulse of junk radiation of about twice the 
amplitude of the earlier run. The new simulation also 
shows smaller extrapolation errors, presumably because 
the new simulation uses larger extraction radii (up to 
r = 350Af, whereas Ref.0] uses a largest extraction ra- 
dius of r = 240 Af). 

We employ four techniques to measure black hole spin: 
Two of these are based on the surface integral for quasi- 
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FIG. 15: Comparison of waveform extrapolation between the 
current simulation of counter-rotating black holes (top panel), 
and the earlier simulation of non-spinning black holes 0, [60| . 
The noise is significantly reduced in the newer simulation, due 
to smaller constraint damping parameters in the wave zone. 



in those regions of the evolution, the black holes can not 
be approximated as isolated Kerr black holes. The be- 
havior of Xsc"(0 9.nd Xgc'^(i) contain information about 
the deformation of the black holes. The final state of 
the simulation is expected to be a single, stationary Kerr 
black hole, for which Xscf (^) Xsc'^(0 should result in 
the correct spin. Indeed, all four spin diagnostics agree 
at very late time to five significant digits (cf. Fig. (5]). 

The accuracy of our simulation places new constraints 
on analytic formulae that predict the final black hole spin 
from the initial spins and masses of a black hole binary. 



Prediction Formula 


Xf 


Mf/M 


Kesden [94] 


0.521153 


0.97039 


Buonanno, Kidder & Lehner [95] 


0.505148 


1.0 


Tichy & Marronetti [96j 


0.548602 


0.962877 


Boyle & Kesden [97] 


0.547562 


0.964034 


Barausse & RezzoUa [98] 


0.546787 


1.0 


Numerical result (this paper) 


0.54781 


0.961109 



TABLE III: Predictions of final black hole spin and mass from 
analytical formulae in the literature, applied to the simulation 
considered here. Refs. [9^, [9^ do not predict the final mass, 
but instead assume zero mass loss. 



local linear momentum, and utilize either simple coor- 
dinate rotation vectors XCoord(i) or approximate Killing 
vectors, XAKv{t); the other two are based on the shape of 
the apparent horizon, and infer the spin from the extrema 
of the scalar curvature (x™(Oj Xsc^i^))- The four spin 
measures agree to better than 1 per cent during the inspi- 
ral. The AKV spin XAKv{t) shows the least variations 
during the simulation, and is the only spin diagnostic 
that results in a monotonically decreasing spin during 
the inspiral, as expected from the effects of tidal fric- 
tion. The other three spin measures (xcooid(Oi X^sc^i^)^ 
Xsc"(0) show various undesircd and physically unreason- 
able behaviors: All three result in increasing spin during 
the inspiral, inconsistent with tidal friction (cf. Fig. d]). 
Xgc"(<) and Xsc^i^)' furthermore show very strong varia- 
tions during the initial transients, just before merger, and 
just after the common horizon forms. This is expected, as 
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